
# Load packages-------
library(tidyverse)
library(dplyr)
library(plm)
library(lmtest)
library(sandwich)
library(glue)
library(haven)
library(ggplot2)
library(ggthemes)
library(broom)

# path & clear
rm(list=ls())
theme_set(theme_few())

# FILL IN PATH TO RUN CODE
dropbox <- "..."
dir <- glue("{dropbox}/_0_bjps_replication")
proc <- glue("{dir}/proc")
graph <- glue("{dir}/out/graph")

# data --------

ISSP_for_analysis <- read_dta(glue("{proc}/ISSP_for_analysis.dta"))

# key variables --------

ISSP_for_analysis$female<-NA
ISSP_for_analysis$female[ISSP_for_analysis$SEX==1]<-0  # Male
ISSP_for_analysis$female[ISSP_for_analysis$SEX==2]<-1  # Female



ISSP_for_analysis$country<-factor(ISSP_for_analysis$country)

ISSP_for_analysis<-ISSP_for_analysis %>%
  mutate(cntry=case_when(country==36~"Australia",
                         country== 40~"Austria",
                         country==152~"Chile",
                         country==191~"Croatia",country==203~"Czech Republic",country==208~"Denmark",
                         country==246~"Finland",country==250~"France",country==276~"Germany",
                         country==352~"Iceland",country==376~"Israel",country==380~"Italy",country==392~"Japan",
                         country==440~"Lithuania",country==554~"New Zealand",country==578~"Norway",
                         country==643~"Russia",country==705~"Slovenia",country==752~"Sweden",
                         country==756~"Switzerland",country==764~"Thailand",
                         country==840~"USA",country==862~"Venezuela"))

ISSP_for_analysis <- ISSP_for_analysis %>%
  mutate(is_OECD = ifelse(country %in% c(36, 40, 152, 203, 208, 246, 250, 276, 352, 380, 392, 554, 
                                         578, 705, 752, 756, 840), 1, 0))


# Analysis --------

lm1_bycountry <- ISSP_for_analysis %>%
  filter(DEGREE>0 & v53>0) %>%
  filter(country %in% c(36,40,203,208,246,250,276,352,376,380,392,554,
                        578,705,752,756,840, 440)) %>%
  group_by(cntry)%>%
  do(tidy(lm(voteright ~SBO+              
               factor(DEGREE)+factor(v53)+factor(female)+AGE,
             weights = WEIGHT, data = .)))


lm_bycountry<-lm1_bycountry%>%
  filter(term=="SBO")

OECD_data <- ISSP_for_analysis %>%
  filter(DEGREE>0 & v53>0) %>%
  mutate(is_OECD = ifelse(country %in% c(36,40,203,208,246,250,276,352,376,380,392,554,
                                         578,705,752,756,840, 440), 1, 0))

lm2_allcountry<-lm(voteright ~  SBO+factor(is_OECD)+
              factor(DEGREE)+factor(v53)+factor(female)+AGE,
            weights = WEIGHT, data =OECD_data)

lm_allcountry<-tidy(lm2_allcountry)     
lm_allcountry<-lm_allcountry%>%
  filter(term=="SBO")%>%
  mutate(cntry = "All Countries")

combined_data <- bind_rows(lm_bycountry, lm_allcountry)
combined_data$cntry <- fct_reorder(combined_data$cntry, desc(combined_data$estimate))
combined_data$cntry <- fct_reorder(combined_data$cntry, combined_data$estimate)


bold_labels <- function(x) {
  ifelse(x %in% c("USA", "All Countries"), "bold", "plain")
}

#  ggplot -----
ggplot <- combined_data %>%
  ggplot(aes(x = cntry, y = estimate, fill = term, colour = term)) +
  scale_fill_manual(values = c("#2F5D62", "#2F5D62")) +
  scale_colour_manual(values = c("#2F5D62", "#2F5D62")) +
  geom_hline(yintercept = 0, linetype = 2, lwd = 0.9, size = 2, colour = "#ACBBBF", alpha = 0.9) +
  geom_errorbar(stat = "identity", alpha = 0.5, position = position_dodge(width = 0.15),
                aes(ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error),
                lwd = 1.5, width = 0) +
  geom_errorbar(stat = "identity", alpha = 0.9, position = position_dodge(width = 0.15),
                aes(ymin = estimate - 1.64 * std.error, ymax = estimate + 1.64 * std.error),
                lwd = 1.3, width = 0) +
  geom_point(stat = "identity", alpha = 0.9, size = 2.5, position = position_dodge(width = 0.15)) +
  coord_flip() +
  #theme_minimal() +
  theme(legend.position = "none",
        axis.text = element_text(size=12),
        axis.text.y = element_text(size = 12, face = bold_labels(levels(combined_data$cntry)))) +
  labs(x = "", y = " ", title = "", colour = "", fill = "", shape = "", group = "")
ggsave(glue("{graph}/issp.png"))
